Introduction to mrgsolve for PBPK

Fall 2022

Lilly, IU, and Metrum Research Group

1 Thank you for inviting us

2 What can you do with mrgsolve?

3 About mrgsolve

4 mrgsolve started as QSP modeling tool

5 Orientation

6 R for Data Science

https://r4ds.had.co.nz/

7 What we will cover today

  1. Three basic workflows
  2. Loading the model into R
  3. Event objects
  4. Data sets
  5. Model specification - code together
  6. Whatever you ask about

Emphasis is on getting you running your own simulations today.

8 A basic simulation with mrgsolve

mod %>% ev(amt = 100, ii = 24, addl = 3) %>% mrgsim() %>% plot()

9 A basic simulation with mrgsolve

mod %>% ev(amt = 100, ii = 24, addl = 3) %>% mrgsim() %>% plot()

10 What’s coming …

11 Why do we use %>% ?

What happens first in this operation?

mean(sqrt(seq(4)))

12 Pipelines

mean(sqrt(seq(4)))
4  %>% seq() %>% sqrt() %>% mean()

Better.

4  %>% 
  seq(.) %>% 
  sqrt(.) %>% 
  mean(., na.rm = TRUE)
mod %>% some_intervention() %>% simulate() %>% post_process()

13 The model object

model %>% ...

14 Take a look: overview

mod


-----------------  source: pk1.cpp  -----------------

  project: /Users/kyleb/ghe...gsolve/models
  shared object: pk1-so-10d6a5b45d9a0 

  time:          start: 0 end: 192 delta: 0.2
                 add: <none>

  compartments:  EV CENT [2]
  parameters:    CL V KA [3]
  captures:      CP [1]
  omega:         0x0 
  sigma:         0x0 

  solver:        atol: 1e-08 rtol: 1e-08 maxsteps: 20k
------------------------------------------------------

15 Take a look: parameters (really important)

param(mod)

 Model parameters (N=3):
 name value . name value
 CL   1     | V    20   
 KA   1     | .    .    

16 Take a look: compartments

init(mod)

 Model initial conditions (N=2):
 name       value . name     value
 CENT (2)   0     | EV (1)   0    

17 Where did mod come from?

mod <- mread("simple", "model")
mod <- mread("<model-name>", "<project-directory>")

18 Read in a model object with caching

First time to read

mod <- mread_cache("simple", here("day1/model"))
. Building simple ... done.

Next time to read

mod <- mread_cache("simple", here("day1/model"))
. Loading model from cache.

19 Model build directory (soloc)

It can be helpful to set

options(mrgsolve.soloc = "my_build_directory")

Equivalent to

mod <- mread("simple", here("day1/model"), soloc = "my_build_directory")

20 Internal model library

mod <- mread("<first-argument>", "<second-argument>")

21 Internal model library

mod <- mread("<first-argument>", "<second-argument>")
mod <- mread("<first-argument>", modlib())
modlib()
. [1] "/Users/kyleb/ghe/courses/pbpk-at-lilly-IU/rlibs/mrgsolve/models"

22 Internal model library

mod <- mread("effect", modlib())
mod
. 
. 
. ----------------  source: effect.cpp  ----------------
. 
.   project: /Users/kyleb/ghe...gsolve/models
.   shared object: effect-so-21f170c05f64 
. 
.   time:          start: 0 end: 36 delta: 0.1
.                  add: <none>
. 
.   compartments:  GUT CENT PERIPH CE [4]
.   parameters:    VC KA K10 K12 K21 E0 EMAX EC50 KEO [9]
.   captures:      EFFECT CP [2]
.   omega:         0x0 
.   sigma:         0x0 
. 
.   solver:        atol: 1e-08 rtol: 1e-08 maxsteps: 20k
. ------------------------------------------------------

23 Try this now

Open R script in Rstudio then …

library(mrgsolve)
mod <- mread("pbpk", modlib())
mod <- modlib("pbpk")

24 Your turn

25 Event objects

e <- ev(amt = 100) 

e
. Events:
.   time amt cmt evid
. 1    0 100   1    1

26 Three ways to invoke

Inline

mod %>% ev(amt = 100) %>% mrgsim()

Object via pipeline

e <- ev(amt = 100)

mod %>% ev(e) %>% mrgsim()

As argument

mod %>% mrgsim(events = e)

27 What to include in ev(...)

28 Interventions and corresponding evid

29 Your Turn

30 Create complex events - 1

What’s going to happen?

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, time = 24, ii = 24, addl = 4)

c(e1, e2)

31 Create complex events - 1

What’s going to happen?

32 Create complex events - 1

Combine two events

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, time = 24, ii = 24, addl = 4)

c(e1, e2)
. Events:
.   time amt cmt evid ii addl
. 1    0 200   1    1  0    0
. 2   24 100   1    1 24    4

33 Create complex events - 2

What’s going to happen?

e1 <- ev(amt = 200, ii = 12, addl = 2) 

e2 <- ev(amt = 100, ii = 24, addl = 4)

seq(e1, e2)

34 Create complex events - 2

What’s going to happen?

35 Create complex events - 2

Put two events in a sequence

e1 <- ev(amt = 200, ii = 12, addl = 2) 

e2 <- ev(amt = 100, ii = 24, addl = 4)

seq(e1, e2)
. Events:
.   time amt ii addl cmt evid
. 1    0 200 12    2   1    1
. 2   36 100 24    4   1    1

36 Create complex events - 3

What is going to happen?

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, ii = 24, addl = 4)

seq(e1, wait = 36, e2)

37 Create complex events - 3

What is going to happen?

38 Create complex events - 3

Wait before starting the next part of the regimen

e1 <- ev(amt = 200) 

e2 <- ev(amt = 100, ii = 24, addl = 4)

seq(e1, wait = 36, e2)
. Events:
.   time amt ii addl cmt evid
. 1    0 200  0    0   1    1
. 2   36 100 24    4   1    1

39 Your Turn

40 Event objects are just data frames

as.data.frame(e1)
.   time amt cmt evid
. 1    0 200   1    1

41 Simulate


model %>% intervention %>% Go!


e <- ev(amt = 100)

mod %>% ev(e) %>% mrgsim()
. Model:  effect 
. Dim:    362 x 8 
. Time:   0 to 36 
. ID:     1 
.     ID time    GUT   CENT PERIPH     CE EFFECT     CP
. 1:   1  0.0   0.00  0.000 0.0000 0.0000  157.0  0.000
. 2:   1  0.0 100.00  0.000 0.0000 0.0000  157.0  0.000
. 3:   1  0.1  91.21  8.443 0.1552 0.2225  155.7  3.460
. 4:   1  0.2  83.19 15.502 0.5817 0.8048  152.8  6.353
. 5:   1  0.3  75.88 21.354 1.2272 1.6382  149.6  8.752
. 6:   1  0.4  69.21 26.156 2.0463 2.6356  146.6 10.719
. 7:   1  0.5  63.13 30.046 3.0001 3.7278  144.1 12.314
. 8:   1  0.6  57.58 33.146 4.0553 4.8607  142.2 13.584
mrgsim(mod, events = e)

42 What would you like to “fix” in this plot?

mod %>% ev(amt = 100) %>% mrgsim() %>% plot(CP~time)

43 Simulation end time

mod %>% ev(amt = 100) %>% mrgsim(end = 48) %>% plot(CP~time)

44 Simulation time step

mod %>% ev(amt = 100) %>% mrgsim(end = 48, delta = 0.1) %>% plot(CP~time)

45 Simulation time grid

stime(mod)
. [1]  0  6 12 18

46 We’re stil working on this setup


model %>% intervention %>% Go! %>% take-a-look


model:

intervention:

47 Working with simulated output


model %>% intervention %>% Go! %>% take-a-look


48 Plot

mod <- mread_cache("effect", modlib())

mod %>% ev(amt = 100) %>% mrgsim() %>% plot()

49 Plot

mod %>% ev(amt = 100) %>% mrgsim() %>% 
  plot(CP + EFFECT~., col = "firebrick")

50 Pipeline to dplyr functions

mod %>% mrgsim() %>% mutate(arm = 1)
. # A tibble: 361 × 9
.       ID  time   GUT  CENT PERIPH    CE EFFECT    CP   arm
.    <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>
.  1     1   0       0     0      0     0    157     0     1
.  2     1   0.1     0     0      0     0    157     0     1
.  3     1   0.2     0     0      0     0    157     0     1
.  4     1   0.3     0     0      0     0    157     0     1
.  5     1   0.4     0     0      0     0    157     0     1
.  6     1   0.5     0     0      0     0    157     0     1
.  7     1   0.6     0     0      0     0    157     0     1
.  8     1   0.7     0     0      0     0    157     0     1
.  9     1   0.8     0     0      0     0    157     0     1
. 10     1   0.9     0     0      0     0    157     0     1
. # … with 351 more rows

51 Controlling output - request

mod %>% mrgsim() %>% head(n = 2)
.   ID time GUT CENT PERIPH CE EFFECT CP
. 1  1  0.0   0    0      0  0    157  0
. 2  1  0.1   0    0      0  0    157  0
mod %>% Req(CP) %>% mrgsim() %>% head(n = 2)
.   ID time CP
. 1  1  0.0  0
. 2  1  0.1  0

52 Controlling output - obsonly

mod %>% ev(amt=1) %>% mrgsim() %>% head(n = 2)
.   ID time GUT CENT PERIPH CE EFFECT CP
. 1  1    0   0    0      0  0    157  0
. 2  1    0   1    0      0  0    157  0
mod %>% ev(amt=1) %>% obsonly() %>% mrgsim() %>% head(n = 2)
.   ID time       GUT       CENT      PERIPH          CE   EFFECT         CP
. 1  1  0.0 0.0000000 0.00000000 0.000000000 0.000000000 157.0000 0.00000000
. 2  1  0.1 0.9121052 0.08443153 0.001551552 0.002224565 156.9866 0.03460309

53 Updating the model object

54 Update

On the fly

mod %>% update(end = 72) %>% mrgsim()

Persistent

mod <- update(mod, end = 72)

mrgsim will call update for you (on the fly)

mod %>% mrgsim(end = 72)

55 What else can I update?

mod %>% update(rtol = 1E-12) %>% ...
mod %>% mrgsim(rtol = 1E-12) %>% ...

56 Updating parameters

mod <- param(mod, CL = 1.5)

57 Add a population element to the simulation

58 One population simulation with mrgsolve

mod %>%
ev(amt = 100, ii = 24, addl = 3) %>%
idata_set(pop) %>%
mrgsim() %>% plot()

59 idata_set takes in individual-level data

head(pop, n = 3)
.   ID   CL   VC     KA KOUT IC50 FOO
. 1  1 1.05 47.8 0.8390 2.45 1.28   4
. 2  2 0.73 30.1 0.0684 2.51 1.84   6
. 3  3 2.82 23.8 0.1180 3.88 2.48   5
length(unique(pop$ID))
. [1] 10

60 pop and mod are connected via parameters

head(pop, n = 3)
.   ID   CL   VC     KA KOUT IC50 FOO
. 1  1 1.05 47.8 0.8390 2.45 1.28   4
. 2  2 0.73 30.1 0.0684 2.51 1.84   6
. 3  3 2.82 23.8 0.1180 3.88 2.48   5
param(mod)
. 
.  Model parameters (N=3):
.  name value . name value
.  CL   1     | V    20   
.  KA   1     | .    .

61 What else can we do with idata?

Batches of simulations or sensitivity analyses

idata <- expand.idata(CL = seq(0.5, 1.5, 0.25))
idata
.   ID   CL
. 1  1 0.50
. 2  2 0.75
. 3  3 1.00
. 4  4 1.25
. 5  5 1.50

62 Note: this is the event + idata_set configuration

mod %>%
  idata_set(idata) %>%
  ev(amt = 100, ii = 24, addl = 2) %>%
  mrgsim(end = 120) %>% plot(log(CP) ~ .)

63 Your turn

64 data_set is the dosing equivalent to idata_set

data <- expand.ev(amt = c(100, 300, 1000), ii = 24, addl = 3)

head(data)
.   ID time  amt ii addl cmt evid
. 1  1    0  100 24    3   1    1
. 2  2    0  300 24    3   1    1
. 3  3    0 1000 24    3   1    1

65 Note: this is data_set configuration

mod %>%
  data_set(data) %>%
  mrgsim(end = 120) %>% plot(log(CP) ~ .)

66 data_set can also carry parameters

data <- expand.ev(
  amt = c(100,300,1000),
  ii = 24, addl = 3,
  CL = seq(0.5,1.5, 0.5)
)

head(data)
.   ID time  amt ii addl cmt evid  CL
. 1  1    0  100 24    3   1    1 0.5
. 2  2    0  300 24    3   1    1 0.5
. 3  3    0 1000 24    3   1    1 0.5
. 4  4    0  100 24    3   1    1 1.0
. 5  5    0  300 24    3   1    1 1.0
. 6  6    0 1000 24    3   1    1 1.0
data <- mutate(data, dose = amt)

mod %>%
  data_set(data) %>%
  mrgsim(carry.out = "dose", end = 120) %>%
  plot(log(CP)~time|factor(dose), group = ID, scales = "same")

67 carry_out or carry.out

Copy from the input data to the output data

mod %>% carry_out(dose) %>% data_set(data, dose==300) %>% mrgsim()
. Model:  pk1 
. Dim:    2886 x 4 
. Time:   0 to 192 
. ID:     3 
.     ID time dose     CP
. 1:   2  0.0  300  0.000
. 2:   2  0.0  300  0.000
. 3:   2  0.2  300  2.712
. 4:   2  0.4  300  4.919
. 5:   2  0.6  300  6.712
. 6:   2  0.8  300  8.167
. 7:   2  1.0  300  9.345
. 8:   2  1.2  300 10.296

68 data_set and ev

head(data, n = 3)
.   ID time  amt ii addl cmt evid  CL dose
. 1  1    0  100 24    3   1    1 0.5  100
. 2  2    0  300 24    3   1    1 0.5  300
. 3  3    0 1000 24    3   1    1 0.5 1000
ev(amt = 100, ii = 24, addl = 3)
. Events:
.   time amt ii addl cmt evid
. 1    0 100 24    3   1    1

69 We have now seen 3 simulation setups

70 Wait a minute …

head(data)
.   ID time  amt ii addl cmt evid  CL dose
. 1  1    0  100 24    3   1    1 0.5  100
. 2  2    0  300 24    3   1    1 0.5  300
. 3  3    0 1000 24    3   1    1 0.5 1000
. 4  4    0  100 24    3   1    1 1.0  100
. 5  5    0  300 24    3   1    1 1.0  300
. 6  6    0 1000 24    3   1    1 1.0 1000

71 Generating data sets

Recall that data sets are just plain old data.frames in R. Feel free to code these up in any way that is convenient for you.

In our experience, the following helper functions cover many (not every) common needs for building the data sets.

72 expand.ev

Like expand.grid

expand.ev(ID = 1:2, amt = c(100,200))
.   ID time amt cmt evid
. 1  1    0 100   1    1
. 2  2    0 100   1    1
. 3  3    0 200   1    1
. 4  4    0 200   1    1

73 ev_rep

e <- ev(amt = 100, ii = 12, addl = 14)
ev_rep(e, ID = seq(5))
.     ID time amt ii addl cmt evid
. 1    1    0 100 12   14   1    1
. 1.1  2    0 100 12   14   1    1
. 1.2  3    0 100 12   14   1    1
. 1.3  4    0 100 12   14   1    1
. 1.4  5    0 100 12   14   1    1
e <- ev(amt = 100, ii = 12, addl = 14, ID = seq(5))

74 ev_rep

e <- seq(ev(amt = 100), wait = 36, ev(amt = 50, ii = 24, addl = 4))
e
. Events:
.   time amt ii addl cmt evid
. 1    0 100  0    0   1    1
. 2   36  50 24    4   1    1
ev_rep(e, ID = seq(2))
.     ID time amt ii addl cmt evid
. 1    1    0 100  0    0   1    1
. 2    1   36  50 24    4   1    1
. 1.1  2    0 100  0    0   1    1
. 2.1  2   36  50 24    4   1    1

75 as_data_set

data <- as_data_set(
  ev(amt = 100, ii = 12, addl = 19, ID = 1:2),
  ev(amt = 200, ii = 24, addl = 9,  ID = 1:3),
  ev(amt = 150, ii = 24, addl = 9,  ID = 1:4)
)
.   ID time cmt evid amt ii addl
. 1  1    0   1    1 100 12   19
. 2  2    0   1    1 100 12   19
. 3  3    0   1    1 200 24    9
. 4  4    0   1    1 200 24    9
. 5  5    0   1    1 200 24    9
. 6  6    0   1    1 150 24    9
. 7  7    0   1    1 150 24    9
. 8  8    0   1    1 150 24    9
. 9  9    0   1    1 150 24    9

76 ev_assign

pop <- expand.idata(GROUP = c(1,2), ID = 1:3)
head(pop)
.   ID GROUP
. 1  1     1
. 2  2     2
. 3  3     1
. 4  4     2
. 5  5     1
. 6  6     2
e1 <- ev(amt = 100, ii = 24, addl = 9)
e2 <- mutate(e1, amt = 200)

data <- ev_assign(
  l = list(e1,e2), 
  idata = pop, 
  evgroup = "GROUP",
  join = TRUE
)

head(data)
.   time amt ii addl cmt evid ID GROUP
. 1    0 100 24    9   1    1  1     1
. 2    0 200 24    9   1    1  2     2
. 3    0 100 24    9   1    1  3     1
. 4    0 200 24    9   1    1  4     2
. 5    0 100 24    9   1    1  5     1
. 6    0 200 24    9   1    1  6     2

77 Your turn

78 Your turn

79 Your turn

80 Model specification

81 Let’s make this model

82 Code up the model

83 Block $PARAM [R] declare and initialize model parameters

$PARAM CL = 1, VC = 20, KA=1.2

84 Block $CMT [txt] declare model compartments

Example $CMT

$CMT GUT CENT

Example $INIT

$INIT GUT = 0, CENT = 10

85 Block $ODE [C++] write differential equations

$ODE
dxdt_GUT  = -KA*GUT;
dxdt_CENT =  KA*GUT - (CL/VC)*CENT;

86 Derive new variables in your model

To get a numeric variable, use double

double x = 5.4;

Other types you might use

bool cured = false;
int i = 2;

If in doubt, use double; it’s what you want most of the time.

87 Block $MAIN [C++] covariate model & initials

For example

$MAIN
double CL = TVCL*pow(WT/70,0.75);
double VC = TVVC*pow(0.86,SEX);
RESP_0 = KIN/KOUT;

In this example

88 C++ expressions and functions

if(a == 2) b = 2;
if(b <= 2) {
  c=3;
} else {
  c=4;
}
double d = pow(base,exponent);
double e = exp(3);
double f = fabs(-4);
double g = sqrt(5);
double h = log(6);
double i = log10(7);
double j = floor(4.2);
double k = ceil(4.2);

Lots of help on the web http://en.cppreference.com/w/cpp/numeric/math/tgamma

89 Preprocessor directives (#define)

$GLOBAL
#define CP (CENT/VC)
#define g 9.8
$ODE
double INH = CP/(EC50+CP);
dxdt_RESP = KIN*(1-INH) - KOUT*RESP;
$CAPTURE CP

90 Introduction to POPULATION simulation

91 Block $OMEGA and $SIGMA [txt]

$OMEGA 0.1 0.2 0.3

92 More $OMEGA and $SIGMA

$OMEGA @block
0.1 0.0947 0.2
$OMEGA @cor
0.1 0.67 0.2

93 Block $OMEGA and $SIGMA

$OMEGA @name PK
0 0 0
$OMEGA @name PD
0 0

Users are encouraged to add labels

$OMEGA @name PK @labels ECL EVC EKA
0 0 0

94 About @ macros

Two forms:

95 Closed-form PK models with $PKMODEL

$PKMODEL cmt = "GUT CENT PERIPH", depot=TRUE
$PARAM CL = 1 , V2 = 30, Q = 8, V3 = 400, KA=1.2

96 More $PKMODEL

$PKMODEL ncmt=1, depot=TRUE
$CMT GUT CENT
$PARAM TVCL = 1 , TVV = 30, KA=1.2
$OMEGA @labels ECL EVC
1 1
$MAIN
double CL = TVCL*exp(ECL);
double V  = TVV *exp(EVC);

97 Bioavailability, dosing lag time, and infusions

98 Functions

add <- function(x,y) {
  z <- x + y
  return(z)
}

add(x = 2, y = 3)
. [1] 5

99 Functions

add <- function(x, y) {
  x + y
}

add(2,3)
. [1] 5

100 Functions

add <- function(x, y=3) {
  x + y
}

add(2)
. [1] 5

101 Sensitivity analysis

102 Local sensitivity analysis

We use the lsa() function from the mrgsim.sa package

library(mrgsim.sa)
?lsa
sens <- lsa(model_object, parameter_names, output_variables)

103 First, look at the result

104 A simulation scenario of interest

mod <- modlib("irm1", end = 72, delta = 0.1)

mod %>%
  ev(amt = 100) %>% 
  mrgsim(end = 72, delta = 0.1) %>% 
  plot()

105 lsa

sens <- lsa(
  mod = ev(mod, amt = 100),
  par = "CL,V2,Q,KA", 
  var = "CP"
)

head(sens, n=3)
. # A tibble: 3 × 5
.    time dv_name dv_value p_name     sens
.   <dbl> <chr>      <dbl> <chr>     <dbl>
. 1   0   CP         0     CL      0      
. 2   0   CP         0     CL      0      
. 3   0.1 CP         0.472 CL     -0.00254

106 The default plot method

plot(sens)

107 Global senstivity analysis

We’ll look at Sobol’s method

108 Global senstivity analysis

109 Global senstivity analysis

110 Sobol